FRidh's blog Discovery consists of seeing what everybody has seen and thinking what nobody has thought -- Albert Szent-Gyorgyi

Articles tagged with work

1. Development of an emission model

   On June the 2nd I will present a paper, titled Determining an Empirical Emission Model for the Auralization of Jet Aircraft, at the Euronoise conference in Maastricht, The Netherlands. My presentation will be in the Auralisation of urban sound session. The conference is just a couple of weeks away from now, and I am still analysing and gathering results. Nevertheless, I thought it would be nice to give a bit of an insight on what I'm working on now and what I will present at Euronoise.

   ## Emission model for auralizations

   Currently I'm developing an emission model for jet aircraft that can be used for auralizations. Existing emission models for noise prediction generally predict sound pressure levels in 1/3-octaves. This works for noise prediction, however, for auralization a finer resolution is needed. Indeed, one needs to be able to model individual tones and in certain cases also modulations.

   At Empa we have in the past conducted measurements for the sonAIR project resulting in a large dataset including audio recordings at multiple sites, cockpit data and flight track data. Colleagues of mine are using this dataset now to develop a next-generation emission model for noise prediction and I'm using this dataset to develop an emission model for auralizations.

   ## Analysis of an event

   Let's now have a look at one specific event and how I analyse such event. I've included the Python code I use for the analysis, so you get a better idea of how I'm working. To give you an idea, there is over 500 GB of audio recordings and several hundreds of MBs on other data. The audio is stored in a single hdf5 file using h5py and all other data in a SQLite database. All data is handled using the amazing Blaze and pandas modules.

   In [1]:

   from sonair.processing import *
   %matplotlib inline
   

   We begin by loading the data belonging to a combination of event (i.e. aircraft passage) and receiver and from a certain start and stop time. start and stop are here seconds relative to a certain event reference time. This reference time is the time at which the aircraft is closest to any of the receivers and turned out to be quite conventient to use.

   In [2]:

   event = '10_004_A320'
   receiver = 'A'
   start = -5.
   stop = +5.
   
   analysis = EventAnalysis(event, receiver, start, stop)
   

   We now have a nice object that gives easy access to all the data. For example, we can request the atmospheric pressure (in mbar) during that event

   In [3]:

   analysis.event.pressure
   

   Out[3]:

   977.39999999999998

   or the coordinates of the receiver (Swiss grid)

   In [4]:

   analysis.receiver.x, analysis.receiver.y, analysis.receiver.z
   

   Out[4]:

   (682692.67099999997, 257054.26800000001, 422.048)

   Obviously, we can also listen to the recording

   In [5]:

   from IPython.display import Audio
   Audio(data=analysis.recording_as_signal, rate=analysis.recording_as_signal.fs)
   

   Out[5]:

   or look at the instantaneous pressure that was measured as a pandas Series

   In [6]:

   analysis.recording.head()
   

   Out[6]:

   2013-09-02 07:25:41.050000       0.618678
   2013-09-02 07:25:41.050022676    0.577433
   2013-09-02 07:25:41.050045352    0.587744
   2013-09-02 07:25:41.050068028    0.659923
   2013-09-02 07:25:41.050090704    0.763036
   Name: signal, dtype: float64

   The above values show you the instantaneous sound pressure as function of time. Perhaps more interesting is to consider the sound pressure level as function of time

   In [7]:

   _ = analysis.recording_as_signal.plot_levels()
   

   

   or the spectrogram, showing the sound pressure not just as function of time, but also as function of frequency.

   In [8]:

   _ = analysis.recording_as_signal.spectrogram(ylim=(0.0, 8000.0))
   

   

   The spectrogram clearly shows the Doppler shift. Also visible, though perhaps not so obvious with the current ranges, is that in the first couple of seconds tones are clearly visible and audible, however, in the latter they're almost gone. This is due to the directivity of the aircraft.

   The emission model needs to be able to predict the emission as function of several parameters, like velocity of the aircraft and/or thrust. What is received by the receiver however also depends on this directivity, i.e., the emission as function of the angles between source and receiver.

   ## Sound propagation

   As the sound waves travel from the source to the receiver several physical phenomena occur. First of all, the aircraft moves with a certain speed relative to the atmosphere, and this causes a Doppler shift. Then, in open spaces sound spreads out, in our case mostly spherically. Furthermore, due to relaxation processes the sound waves are attenuated. The receiver could move as well resulting in an additional Doppler shift. Luckily, we were using fixed microphones, so this is a problem we won't have to deal with.

   However, an unwanted effect that can occur for which correcting would be very difficult are contributions due to reflections. And unfortunately we had ground reflections. To make matters worse, atmospheric turbulence adds a bit of randomness to the propagation. Joy!

   ## Calculating back to the source

   In order to determine the actual emission of the aircraft we now have to undo these effects. This, however, is not that easy. I won't go into details regarding assumptions I made, nor the errors I would make, but instead just briefly explain what I did.

   Using meteorological information I could estimate the speed of sound and attenuation in the atmosphere. Together with positional information of the aircraft and receiver I could undo the time lag and thus Doppler shift (basically by doing the inverse of what I posted earlier on the Doppler shift) as well as attenuation due to air absorption. These operations are all applied directly to the audio recording. Let's listen to the result!

   In [9]:

   Audio(data=analysis.reverted_as_signal, rate=analysis.reverted_as_signal.fs)
   

   Out[9]:

   This sounds quite different from the original recording! The Doppler shift is gone, as you can also see in the spectrogram

   In [10]:

   _ = analysis.reverted_as_signal.spectrogram(ylim=(0.0, 8000.0))
   

   

   and the level also doesn't change that much anymore over time either.

   In [11]:

   _ = analysis.reverted_as_signal.plot_levels()
   

   

   Unfortunately there was a click audible at the start. Fixing that is on my to-do list.

   ## Investiging the emission spectrum

   Now that we've calculated back to the source it is time to investigate the emission. The power spectral density shows many peaks, and there seems to be a pattern in the recurrence. It looks as if there is some constant spacing between the tones.

   In [12]:

   fig = analysis.reverted_as_signal.plot_power_spectrum(xscale='linear', xlim=(500, 3500))
   

   

   Maybe we can check if this is the case and determine this spacing? Luckily, we can. Using cepstral analysis, or more specifically, using the complex cepstrum, we can determine the spacing between these recurrences. I won't show any further how I use the complex cepstrum but instead go straight to the result. The spacing corresponds to the frequency at which the engine's shaft rotates. This frequency is the $N_1$ frequency and is determined to be

   In [13]:

   analysis.N1
   

   Out[13]:

   73.13432836293677

   The tones that you see and hear are multiples of this $N_1$ frequency and are called Buzz-Saw tones. However, one very strong peak is visible at

   In [14]:

   blades = 36
   analysis.N1 * blades
   

   Out[14]:

   2632.8358210657238

   This is the blade passing frequency, which is obtained by multiplying $N_1$ with the amount of blades of the fan. Multiples of the blade passing frequency are also present and audible. The figure below shows again the power spectral density along with black and red lines representing respectively Buzz-Saw tones and multiples of the blade passing frequency.

   In [15]:

   fig = analysis.reverted_as_signal.plot_power_spectrum(xscale='linear', xlim=(1500, 6000))
   for harmonic in analysis.N1_harmonics:
       fig.axes[0].axvline(harmonic, color='black')
   for harmonic in np.arange(1, 3)*analysis.N1*blades:
       fig.axes[0].axvline(harmonic, color='red')      
   

   

   ## Extracting features

   The next step is to estimate the power of the tones and the noise for all events and receiver combinations and using those estimations along with estimations of velocity and thrust to develop an emission model. As mentioned before, this means taking into account the directivity and thus the orientation from source to receiver (or vice versa). The figure below shows the components of the unit vector pointing from receiver to source.

   In [16]:

   _ = analysis.orientation_receiver_to_source.ix[::10000].plot()
   

   

   Right now I'm writing the code needed to automatically extract the features. Hopefully, in a couple of days I can begin analysing the extracted features. I'm looking forward to that!

   ## Then what?

   After I've extracted the features and developed the initial model I plan on doing some preliminary listening tests in our lab at Empa this summer. Depending on the results of those tests, I plan further listening tests in autumn at SP's listening room, which features a 3rd order Ambisonics system. And after that? Well, then I think it's time to write my thesis!

2. SONORUS meeting in Naples

   Two weeks ago there was a SONORUS meeting and school in Naples and Sorrento, Italy. The topic of the school was 'Computational Soundscape Analysis'. The next meeting will be in October in Gothenburg, Sweden.

3. Paper at Euronoise

   The timetable for Euronoise 2015 is available. I will present the paper Determining an Empirical Emission Model for the Auralization of Jet Aircraft on Tuesday the second at nine in the morning.

   Update

   Seems like my presentation is moved to 10:20.

4. Slides on Zenodo

   As posted earlier I've been uploading to Zenodo almost all slides, posters and proceedings that I made during my project. This includes the slides (as well as two auralizations) of the presentation I gave at the 168th ASA meeting in Indianapolis last year.

5. Auralization: Doppler effect

   The goal of the project that I'm working on is to develop a tool to synthesize how it sounds when an aircraft flies over an urban environment. One part of the project deals with the propagation model, which also includes the Doppler shift. Let's have a look now at how we can include the Doppler shift in auralizations.

   Doppler shift equation¶

   The following equation is likely known to you. $$ f = \frac{c + v_r}{c + v_s} f_0 $$ You put in the emission frequency $f_0$, the speed of sound $c$, and the velocity of respectively receiver and source $v_r$ and $v_s$. What you get out is the Doppler shifted frequency.

   In the time-domain the Doppler shift is implicitly included when considering the propagation delay between source and receiver. It takes $d = s/c$ seconds to travel the source-receiver distance $s$ at the speed of sound $c$. When the distance changes over time, the emission is Doppler shifted.

   Applying the Doppler shift to a tone¶

   As an example we will now apply the Doppler shift to a pure tone. Let's start of with some imports.

   In [1]:

   import numpy as np
   %matplotlib inline
   import matplotlib.pyplot as plt
   from acoustics import Signal
   from IPython.display import Audio
   import warnings
   warnings.filterwarnings('ignore')
   

   As you can see I'm importing the Signal class from the acoustics package. This package is mainly a collection of tools that could be useful to acousticians, and I started this package because I need many of these tools for my research, and I'm sure others can use it as well. The code base is still a bit messy though, and things do break occasionally.

   We start with generating a pure tone.

   In [2]:

   duration = 10.0
   fs = 44100.0
   frequency = 5000.0
   
   samples = int(fs*duration)
   times = np.arange(samples) / fs
   
   signal = Signal(np.sin(2.0*np.pi*frequency*times), fs=fs)
   

   As I wrote before the Signal class can be found in the acoustics package. This class is a subclass of the np.ndarray and requires besides the array also a sample frequency. The class has some methods for common operations, like calculating different types of spectra and plotting them, as well as calculating other acoustic quantities like the equivalent sound pressure level.

   In [3]:

   Audio(signal, rate=signal.fs)
   

   Out[3]:

   In [4]:

   fig = signal.spectrogram(ylim=(0.0, 10000.0))
   

   

   In order to create a Doppler shift we need a moving source or receiver. In this case, let's consider a moving source and a non-moving receiver. The source will move straight over the receiver.

   In [5]:

   velocity = 70.0
   
   x = times * velocity
   x -= x.max()/2.0
   y = np.zeros(samples)
   z = 100.0 * np.ones(samples)
   
   position_source = np.vstack((x,y,z)).T
   position_receiver = np.zeros(3)
   

   In [6]:

   distance = np.linalg.norm((position_source - position_receiver), axis=-1)
   

   The exact speed of sound in the atmosphere depends on conditions of the atmosphere, like the temperature, pressure, humidity. Furthermore, the speed of sound can vary spatially and temporally due to e.g. wind or temperature gradients, or turbulence. Let's assume for now a constant speed of sound of 343.0 m/s.

   In [7]:

   soundspeed = 343.0
   

   The propagation delay from source to receiver is given simply by

   In [8]:

   delay = distance / soundspeed
   

   Let's have a look at the delay as function of time.

   In [9]:

   fig = plt.figure()
   ax = fig.add_subplot(111)
   ax.plot(times, delay)
   ax.set_xlabel('$t$ in s')
   ax.set_ylabel('delay in s')
   

   Out[9]:

   <matplotlib.text.Text at 0x7f77211c85f8>

   

   Now that we've calculated the delay we need to apply this delay to our signal of interest, which is a matter of shifting the samples. However, since we have a discrete signal, and the delays are not integer multiples of the sample time, an interpolation scheme is needed.

   Linear interpolation¶

   The easiest and likely fastest method is a linear interpolation.

   In [10]:

   def interpolation_linear(signal, times, fs):
       """Linear interpolation of `signal` at `times`.
       
       :param signal: Signal.
       :param times: Times to sample at.
       :param fs: Sample frequency.
       
       """
       k_r = np.arange(0, len(signal), 1)        # Create vector of indices
       k = k_r - times * fs                      # Create vector of warped indices 
       kf = np.floor(k).astype(int)              # Floor the warped indices. Convert to integers so we can use them as indices.
       R = ( (1.0-k+kf) * signal[kf] + (k-kf) * signal[kf+1] ) * (kf >= 0) #+ 0.0 * (kf<0)
       return R
   

   Let's apply the delay and listen.

   In [11]:

   delayed_linear = Signal( interpolation_linear(signal, delay, signal.fs), signal.fs)
   

   In [12]:

   Audio(delayed_linear, rate=delayed_linear.fs)
   

   Out[12]:

   You can hear the Doppler shift clearly but...there are artifacts! The artifacts are also clearly visible in the spectrogram shown below. The red line shows the Doppler shifted frequency and the other lines that are visible are the artifacts. The artifacts are at least 30 dB less strong than the tone of interest, but we can still clearly hear them since there is no signal besides the tone.

   In [13]:

   fig = delayed_linear.spectrogram(ylim=(0.0, 10000.0))
   

   

   In most auralizations I create these artifacts are masked. However, there are ways to reduce them. One way would be to upsample the signal before applying this resampling. Another way would be to consider another interpolation scheme.

   Lanczos interpolation¶

   Theoretically, the optimal reconstruction filter is a sinc filter. In practice only approximations of this filter can be used, and these approximations are generally achieved by windowing and truncating the sinc function. One of these approximations is the Lanczos filter or kernel, which is the sinc function windowed by another sinc function. For more information on Lanczos interpolation I can recommend the excellent Wikpedia article on this topic.

   The following code shows an implementation of Lanczos interpolation in one dimension. In pure Python this would be too slow so I'm using here the Numba just-in-time compiler.

   In [14]:

   import numba
   import math
   
   @numba.jit(nogil=True)
   def sinc(x):
       if x == 0:
           return 1.0
       else:
           return math.sin(x*math.pi) / (x*math.pi)
   
   @numba.jit(nogil=True)
   def _lanczos_window(x, a):
       if -a < x and x < a:
           return sinc(x) * sinc(x/a)
       else:
           return 0.0
   
   @numba.jit(nogil=True)    
   def _lanczos_resample(signal, samples, output, a):
       """Sample signal at float samples.
       """
       for index, x in enumerate(samples):
           if x >= 0.0 and x < len(signal):
               for i in range(math.floor(x)-a+1, math.floor(x+a)):
                   if i >= 0 and i < len(signal):
                       output[index] += signal[i] * _lanczos_window(x-i, a)
       return output
   
   
   def interpolation_lanczos(signal, times, fs, a=10):
       """Lanczos interpolation of `signal` at `times`.   
           
       :param signal: Signal.
       :param times: Times to sample at.
       :param fs: Sample frequency.
       :param kernelsize: Size of Lanczos kernel :math:`a`.
       
       """
       samples = -times * fs + np.arange(len(signal))
       return _lanczos_resample(signal, samples, np.zeros_like(signal), a)
   

   Let's apply again the delay but now using this filter. Note that I use the default kernelsize, which is quite arbitrary.

   In [15]:

   delayed_lanczos = Signal(interpolation_lanczos(signal, delay, signal.fs), signal.fs)
   

   In [16]:

   Audio(delayed_lanczos, rate=delayed_lanczos.fs)
   

   Out[16]:

   The artifacts that were present with linear interpolation are mostly gone now, however, when the change in delay is large, other artifacts are audible.

   In [17]:

   fig = delayed_lanczos.spectrogram(ylim=(0.0, 10000.0))
   

   

   ## Conclusions

   The Doppler effect is an interesting physical phenomenon, one that we hear perhaps daily. I've shown how you can synthesize the Doppler shift and considered two interpolation schemes. Which method is preferable depends on your use case. For higher frequencies Lanczos interpolation results in fewer artifacts than linear interpolation. However, computationally Lanczos interpolation is much more demanding:

   In [18]:

   %timeit interpolation_linear(signal, delay, signal.fs)
   %timeit interpolation_lanczos(signal, delay, signal.fs)
   

   10 loops, best of 3: 26.7 ms per loop
   1 loops, best of 3: 597 ms per loop
   

6. Research visit to Chalmers

   After the ASA meeting in Indianapolis I went to Gothenburg again for a research visit. As I wrote in my previous post it was strange to be there again. Last spring I visited Gothenburg for just a couple of days, but now, staying for seven weeks again is quite different. You start noticing quite a lot of things have changed. Many people I knew there have moved away and I also noticed a couple of restaurants I visited have closed since. Despite those things it was nice meeting old friends again.

   So, as I wrote in the title I went there on a research visit. The main reason for my visit was to attend some courses at Chalmers but to also work on my research. Initially the idea was to investigate further the model we're developing to deal with atmospheric turbulence in aircraft auralisations. However, in the end the courses took more time then I expected, so I didn't really get to work on the turbulence model anymore. However, I did manage to implement one nice part of the auralisation tool, namely support for Ambisonics.

   Besides Chalmers I also visited SP in Borås. SP was a familar place for me as well since I did an internship there in 2010. Earlier I had made recordings with a SoundField microphone, however, at Empa we do not have an Ambisonics setup. SP has a listening room with a third order Ambisonics system, giving me the possibility to listen to the first order recording as well as to higher order auralisations.

   One visit I met up with Peter Lundén who installed the system. After some problem solving we managed to get both the recordings and auralisations working. Unfortunately there was one problem with the auralisations, somehow halfway the aircraft makes a U-turn. While I haven't checked it any further, I think it is due to the inclusion of the Condon-Shortley phase in the function I used to calculate the spherical harmonics.

7. Forum Acusticum 2014, Krakow

   Last week the conference Forum Acusticum 2014 took place in the city of Krakow, Poland. From Monday till Friday there were presentations as well as meetings of technical committees. The presentations covered a wide range of topics like structural acoustics, human response and auralisation.

   A special session was organised on Urban Sound Planning as well, in which many of my project members presented their current work. I'm glad a lot of people attended the session.

   I gave a presentation with the title Modelling sound propagation in the presence of turbulence for the auralisation of aircraft noise in a session on Environmental Acoustics. This session was Friday morning early, which was a bit unfortunate, considering there was a nice party with the Young Acousticians Network the evening before. I was surprised there were still so many people that morning!

   On Thursday I joined a discussion on benchmarking for computational acoustics. When implementing a model you would like to be able to validate this. This can sometimes be difficult due to lack of available information or benchmarks. Therefore, the idea is now to put online some kind of repository where people can put their benchmarks. Considering this is a topic of interest in other fields of science I am really curious how they deal with this.

   The city of Krakow surprised me a lot. It's quite a green city and there seems to be a good and relaxed mood. There's plenty of pubs and restaurants, especially around the main square which is huge.

   Anyway, that's my summary on FA2014!

8. Summer in Norway

   Many that know me are aware of the fact that I am spending this summer in Oslo, Norway. In May I was offered a position for the summer at the Noise & Vibration section at Det Norske Veritas in Høvik. Thus far I am quite content with the work. I've been traveling a lot across Norway for noise and vibration measurements (DNV Comfort Class) and I've started working on noise and vibration predictions.

   Oslo seems to be a nice city as well. Especially the location where I live; everything is basically around the corner. Commuting time is acceptable and nature can be reached in about half an hour. Unfortunately I haven't really been able to see much of the nature outside Oslo due to the bad weather and a bad fever last week.

   The coming week I am having off from work, and what way to spend it better than to go hiking! Monday I will take the bus/train to Liomseter from where I will hike for five days along the final part of Jotunheimstien. Last week I bought all the gear (ouch, Gore-Tex is expensive!!) I think I will need for this and future trips. What future trips, you might ask. Well, at least end of next month I intend to walk the Kungsleden in Sweden.

   Anyhow, summer thus far has been quite nice with visiting The Netherlands again, as well as attending Hultsfredsfestivalen in Sweden. And with the forecast of the hiking trips it only seems to become better :-)

   To end this post I have here a nice recording of the Norwegian band Katzenjammer, who were playing at Hultsfredsfestivalen as well, though this is a recording from elsewhere.

Page 1 / 2 » ⇉